function dif = HoaDiffuseness(hoaSig,hoaFmt)
% HOADIFFUSENESS     Diffuseness of a HOA sound scene
%
% dif = HoaDiffuseness(hoaFmt,hoaSig)
%
% Calculate the diffuseness index of a HOA sound scene by looking at the
% eigenvalues of the HOA signal covariance matrix.
%
% Inputs: - hoaSig is a N x M matrix, where N is the number of time samples
%           (or frequency bins) and M is the number of HOA signals.
%         - (optional) hoaFmt defines the format of the HOA signals
%           (typically generated by the GenerateHoaFmt function). If no
%           format is passed it is assumed that the signals are normalised
%           using the 'N2D' or 'N3D' conventions.
%                
% Output: - dif is the diffuseness index for the input set of signals, 
%           comprised between 0 and 1.

% Convert HOA signals to N3D/N2D format if necessary
if nargin > 1
    if ~strcmp(hoaFmt.conv,'N3D') || ~strcmp(hoaFmt.conv,'N2D')
        % Create target hoaFmt structure
        hoaFmtNrm = hoaFmt ;
        if hoaFmt.res3d == 0
            hoaFmtNrm.conv = 'N2D' ;
        else
            hoaFmtNrm.conv = 'N3D' ;
        end
        % Generate hoaFmt conversion matrix
        CnvMat = sparse(Hoa2HoaConversionMatrix(hoaFmt,hoaFmtNrm)) ;
        % Apply normalisation to the HOA signals
        hoaSig = hoaSig * CnvMat.' ;
    end
end

% Number of spherical harmonics
nmbHrm = size(hoaSig,2) ;

% HOA signal covariance matrix
CovMat = cov(hoaSig) ;

% Covariance matrix eigenvalue
eigVal = eig(CovMat) ;

% "variance" of the eigenvalues
eigVar = sum(abs(eigVal-mean(eigVal)))/sum(eigVal) ;

% Diffuseness
dif = abs( 1 - eigVar * nmbHrm/(nmbHrm-1)/2 ) ;
